C Boundary layer flow

C Boundary layer flow past a flat plate
CNumerical method: finite difference

      PARAMETER (NN=1000)
      COMMON/MSH1/NEX,NEY,NE
      COMMON/MSH2/NX,NY
      DIMENSION Y(NN)
      DIMENSION UO(NN),VO(NN)
      DIMENSION U(NN),V(NN)
      DIMENSION F(NN),A(NN,NN)

c  NX: TOTAL NUMBER OF NODES IN X-DIRECTION
c  NX: TOTAL NUMBER OF NODES IN Y-DIRECTION


c  X: COORDONATE
c  Y: COORDONATE


c  ENTER DATA

      WRITE(6,*) 'ENTER STEP SIZE IN X-DIRECTION'
      READ(5,*) DX
      WRITE(6,*) 'ENTER STEP SIZE IN Y-DIRECTION'
      READ(5,*) DY
      WRITE(6,*) 'ENTER NX '
      READ(5,*) NX
      WRITE(6,*) 'ENTER NY'
      READ(5,*) NY
      WRITE(6,*) 'ENTER RHO'
      READ(5,*) RHO
      WRITE(6,*) 'ENTER AMU'
      READ(5,*) AMU
      WRITE(6,*) 'ENTER U FREESTREAM'
      READ(5,*) U8

c    CREAT MESH

      DO J=1,NY
      Y(J)=(J-1)*DY
      END DO

c    SET BPOUNDARY CONDITIONS

      PI =ACOS(-1.)
      NH=2*NY/3
      NP=2*NY

      DO J=1,NH
      VO(J)=0.0
      UO(J)=U8*SIN(PI/2.0*(J-1)/(NH-1))
      END DO
      DO  J=NH+1,NY
      VO(J)=0.0
      UO(J)=U8
      END DO

      DO I=1,NY
      WRITE(6,*) I,Y(I),UO(I),VO(I)
      END DO

c     SET INITIAL GUESS
      DO  J=1,NY
      V(J)=VO(J)
      U(J)=UO(J)
      END DO

      ISTEP=0
      ERRMAX=1.E-5
10    ISTEP= ISTEP+1

c     START ITERATION PROCEDURE

      ITER=0
20    ITER=ITER+1

      DO I=1,NP
      DO J=1,NP
      A(I,J)=0.0
      END DO
      END DO

      DO I=1,NY
      IF(I.EQ.1) THEN
      F(1)=0.0
      F(2)=0.0
      A(1,1) =1.0
      A(2,2) =1.0
      END IF

      IF (I.EQ.NY) THEN
      F(NP-1)=(V(I)-V(I-1))/DY
      F(NP)=0.0

      A(NP-1,NP-3)=-1.0/DY
      A(NP-1,NP-1)=1./DY
      A(NP,NP)=1.0

      END IF

      IF ((I.LT.NY).AND .(I.GT.1))THEN
      F(2*I-1)=(U(I)-UO(I))/DX+(V(I+1)-V(I-1))/(2.0*DY)
      F(2*I)=U(I)*(U(I)-UO(I))/DX+V(I)*(U(I+1)-U(I-1))/(2.0*DY)
      F(2*I)=F(2*I)-AMU/RHO*(U(I+1)-2.0*U(I)+U(I-1))/DY**2

      A(2*I-1,2*I-3)= -0.5/DY
      A(2*I-1,2*I-2)= 0.0
      A(2*I-1,2*I-1)= 0.0

      A(2*I-1,2*I)= 1./DX
      A(2*I-1,2*I+1)= 0.5/DY
      A(2*I-1,2*I+2)= 0.0

      A(2*I,2*I-3)= 0.0
      A(2*I,2*I-2)= -AMU/RHO/DY**2-0.5*V(I)/DY
      A(2*I,2*I-1)= 0.5*(U(I+1)-U(I-1))/DY
      A(2*I,2*I)= (U(I)-UO(I))/DX+U(I)/DX+2.*AMU/RHO/DY**2
      A(2*I,2*I+1)= 0.0
      A(2*I,2*I+2)=-AMU/RHO/DY**2+0.5*V(I)/DY

      END IF

      END DO
C    SOLUTION FOR INCREMENTS

      CALL GAUSS(A,F,NP,NN)

C     UPDATE VARIABLES

      ERRORS=0.0
      DO I=1,NY
      DV=F(2*I-1)
      IF(ABS(DV).GT.ERROR)ERROR=ABS(DV)
      IF(I.EQ.1)DV=0.0
      V(I)=V(I)-DV
      DU=F(2*I)
      IF(ABS(DU).GT.ERROR)ERROR=ABS(DU)
      IF((I.EQ.1).OR.(I.EQ.NY)) DU=0.0
      U(I)=U(I)-DU
      END DO
      WRITE(6,*)'ITER=',ITER,' ERROR=',ERROR
      IF(ERROR.GT.ERRMAX)GOTO 20

C    PRINT SOLUTION, ADVANCE SOLUTION DOWNSTREAM
      DO I=NY,1,-1
      WRITE(6,*),I,Y(I),U(I),V(I)
      UO(I)=U(I)
      VO(I)=V(I)
      END DO
      WRITE(6,*)' ISTEP= ', ISTEP
      IF(ISTEP.LT.NX) GOTO 10
      STOP
      END

      SUBROUTINE GAUSS(A,B,N,NP)
      PARAMETER (NMAX=5000)
      DIMENSION A(NP,NP),B(NP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)

      DO 11 J=1,N
      IPIV(J)=0
11    CONTINUE
      DO 22 I=1,N
      BIG=0
      DO 13 J=1,N
      IF (IPIV(J).EQ.1) THEN
      DO 12 K=1,N
      IF (IPIV(K).EQ.0) THEN
      IF (ABS(A(J,K)).GE.BIG) THEN
      BIG=ABS(A(J,K))
      IROW=J
      ICOL=K
      ENDIF
      ELSE IF(IPIV(K).GT.1) THEN
      PAUSE'Singular matrix'
      ENDIF
12    CONTINUE
      END IF
13    CONTINUE
      IPIV(ICOL)= IPIV(ICOL)+1
      IF (IROW.NE.ICOL) THEN
      DO 14 L=1,N
      DUM=A(IROW,L)
      A(IROW,L)=A(ICOL,L)
      A(ICOL,L)=DUM
14    CONTINUE

      DUM=B(IROW)
      B(IROW)=B(ICOL)
      B(ICOL)=DUM

      ENDIF
      INDXR(I)=IROW
      INDXC(I)=ICOL
      IF (A(ICOL,ICOL).EQ.0. ) PAUSE'SINGULAR MATRIX'
      PIVINV=1./A(ICOL,ICOL)
      A(ICOL,ICOL)=1.
      DO 16  L=1,N
      A(ICOL,L)=A(ICOL,L)*PIVINV
16    CONTINUE

      B(ICOL)=B(ICOL)*PIVINV

      DO 21 LL=1,N
      IF (LL.NE.ICOL) THEN
      DUM=A(LL,ICOL)
      A(LL,ICOL)=0.
      DO 18 L=1,N
      A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18    CONTINUE

      B(LL)=B(LL)-B(ICOL)*DUM

      ENDIF
21    CONTINUE
22    CONTINUE

      DO 24 L=N,1,-1
      IF (INDXR(L).NE.INDXC(L)) THEN
      DO 23 K=1,N
      DUM=A(K,INDXR(L))
      A(K,INDXR(L))=A(K,INDXC(L))
      A(K,INDXC(L))=DUM
23    CONTINUE
      ENDIF
24    CONTINUE
      RETURN
      END